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We show that a laser beam which propagates through an optical medium with Kerr (focusing) 
and higher order (defocusing) nonlinearities displays pressure and surface-tension properties yielding 
capillarity and dripping effects totally analogous to usual liquid droplets. The system is reinterpreted 
in terms of a thermodynamic grand potential, allowing for the computation of the pressure and 
surface tension beyond the usual hydrodynamical approach based on Madelung transformation and 
the analogy with the Euler equation. We then show both analytically and numerically that the 
stationary soliton states of such a light system satisfy the Young-Laplace equation, and that the 
dynamical evolution through a capillary is described by the same law that governs the growth of 
droplets in an ordinary liquid system. 



PACS numbers: 03.75.Lm, 42.65.Jx, 42.65.Tg 

Introduction.- Since the pioneering paper of Piekara 
in the 70's|l), many works have highlighted the inter- 
esting properties of laser beams whose propagation is 
described by the so-called cubic-quintic (CQ) nonlinear 
Schrodinger equation (NLSE) with competing nonlin- 
earities. Cavitation, superfluidity and coalescence have 
been investigated [2, 3] in the context of liquid He, where 
the model is a simple approach which does not take 
into account nonlocal interactions. Stable optical vor- 
tex solitons and the existence of top-flat states have 
also been reported in optical materials with CQ optical 
susceptibility [4]. Recent experiments about filamenta- 
tion of high-power laser pulses have shown that the CQ 
regime could be achievable in CS2 [5] as well in some 
chalcogenide glasses [6j. Recently, it has been also sug- 
gested that atomic coherence may be used to induce a 
giant CQ-like refractive index in Rb gas [7|]. 

On the other hand, this model has been shown to dis- 
play surface properties in numerical simulations of soliton 
collisions [8], that have been considered as a trace of a 
liquid state of light [7|, |9|. Here, we will provide the first 
analytical, quantitative demonstration of the liquid be- 
havior of the system in (2 + 1) dimensions, both in its 
stationary soliton solutions and in the dynamical evo- 
lution when a light bump is forced to pass through a 
wave-guide simulating a capillary. In particular, we will 
provide the first consistent computation of the pressure 
and surface tension of the light bubbles in (2 + 1) di- 
mensions, and show both analytically and numerically 
that they satisfy the Young-Laplace (Y-L) equation that 
gives the equilibrium of usual droplets. Subsequently, we 
will demonstrate that the system dripping properties are 
governed by the same generalized Y-L that applies to an 
ordinary liquid. These results show the deep connection 
between the nonlinear dynamics of laser beams and co- 
herent liquids at zero temperature [lOj . 

Thermodynamic model- We will consider the paraxial 
propagation through an ideal CQ medium of a linearly 
polarized laser beam, being its complex amplitude distri- 



bution \I> described by a nonlinear Schrodinger equation 
of the adimensional form: 
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where z is the propagation distance multiplied by 27r/A, 
being A the wavelength of the continuous light beam; 
X7 2 ± is the transverse Laplace operator in terms of x,y, 
the spatial variables multiplied by 2i\yJ2n^j\ being no 
the linear refractive index of the medium; 7 and S are 
proportional to the (opposite) an d optical sus- 
ceptibilities, respectively. 

It is well-known that stationary version of Eq.([Tj) 
admits localized soliton-like solutions [1] of the form 
*&a(x, V: z) = y)e~ ZfJjZ , being fi the propagation con- 
stant. In particular, it has been shown numerically that 
high power solitons feature top-flat profiles [4]. These 
modes can only be calculated numerically and coex- 
ist with plane waves solutions of constant amplitude 
^fA(x,y,z) = Ae~ ZfJjZ , which lead by substitution in 
Eq.([T]) to = 5\A\ 4 — j\A\ 2 . The existence domain for 
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solitons is fi G (/ioo,0), where fi^ = — 

As discussed in [11], the stationary solutions of Eq. (pp) 
can be derived from a variational principle J^pr = from 
the Landau's grand potential Vt = H — fiN, where H is 
the Hamiltonian, N — J \ ^\ 2 dS the particle number (in 
our system the photon flux through the transverse section 
S) and /i the chemical potential which is thus identified 
in our case as the propagation constant we have defined 
above. The resulting expression is: 
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dS, (2) 



In three dimensional systems, the partial derivative 
(dVt/dV)^ T of Vt with respect to the volume V at con- 
stant chemical potential /i and temperature T would give 
minus the pressure. Assuming a completely coherent 



2 



two-dimensional model with T = OK, it is then natu- 
ral to use the derivative of Q with respect to the area <S, 
which yields: 
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FIG. 1: [Color online] Plot of radius vs. propagation constant 
for different localized stationary solutions of Eq.(pQ). Solid 
(dashed) line corresponds to numerical (analytical) calcula- 
tions. The left (right) inner picture, shows the field modu- 
lus profile (pressure distribution) of a stationary beam with 
/V/ioo — 0.98. Notice that close to the origin, the homo- 
geneous pressure is positive whereas in the inhomogeneous 
region becomes negative. 

In Fig [1] we plot an example of a "flat-top" stationary 
state, corresponding to p/ Poo = 0.98, together with its 
pressure distribution. For all this type of solutions, in the 
wide flat region the value of the pressure is positive and 
equal to the central value p{0) = p c . Close to the origin 
it is straightforwardly obtained for the beam amplitude: 



|A(0,0)| 2 = 
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which yields for the pressure at the center of the beam 
(p c ) the following analytical expression: 
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Notice that p c vanishes at the limits of the existence 
domain: p = (trivial zero-amplitude solution) and 
p = /ioo (infinite "flat-top" with {A^ 2 = 0.757/(5). This 
result can be considered as a trace of the existence of two 



possible vacuum states in the system [11]. 

Limitations of the hydrodynamical analogy.- Our 
present purpose is to use the previous formalism in or- 
der to study the physical properties of the system. First, 
let us try to apply the commonly used hydrodynamical 
analogy, based on introducing the so-called Madelung 
Transformation [l2| ip = p 1 / 2 e^, where p is a positive 
definite real function and (j) is a real phase, eventually 
depending on the variables x, y and z. After substi- 
tuting in Eq. (pp), and separating both the real and the 



imaginary part, we get two equations. The first one is a 
continuity equation, that is used to establish an analogy 
with hydrodynamics, and to argue that can be iden- 
tified with a current, i.e. a velocity field v. The second 
equation is 
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Here, the usual approach [2] is to assume that the term 
2p i/2 V^p 1 / 2 can be neglected. In this case, Eq. (|6|) is 
identical to the known Euler equation of hydrodynamics 
provided that \/±(jp — Sp 2 ) = — This formula was 
used e.g. in Ref. [2] in order to derive an expression 
for the pressure of the homogenous phase, which coin- 
cides with our result of Eq. ([5]) for the top flat region 
of the solitons. One could hope that this analogy could 
be generalized also to the non-homogenous zones, such 
as that corresponding to values of the radial coordinate 
r around the radius R of the beam. However, we will 
see that this is not the case. For a (non-rotational) top 
flat soliton the phase term is simply equal to (j) = —pz, 
therefore v = V±(j) = 0. Note that p is the propagation 
constant which is just a constant number for the given 
soliton. Therefore, by substituting in Eq. (|6]), we get 
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In other words, in this case the term that is usually 
neglected is exactly equal and opposite to the part that is 
used to compute the pressure. Therefore, neglecting such 
a term would correspond to a 100% error. In fact, for 
any soliton solution, having v = 0, Euler equation would 
unavoidably imply a constant pressure, which cannot be 
the case in any region where the spatial variation of the 
density is important. Nevertheless, in spite of this failure 
of the ordinary hydrodynamical approach, we will see 
that the 'pressure' distribution in the non- homogeneous 
region can still be given a deep physical interpretation. 

Equilibrium and surface tension.- Fig. 1 suggests that 
the Q potential, as given by the spatial integral of minus 
the pressure, can be expressed as the sum of two contri- 
butions: one from the top-flat region, which is —7rR 2 p c 
with a good accuracy; and another from the region near 
the border, where the field is spatially-dependent, which 
is — 2tt rp(r)dr. Thus, we get the following analyti- 
cal approximation for the Q potential of the "flat-top" 
solutions: 



-7rR 2 p c + 2naR, 



(8) 



where we have defined a parameter a = — rp(r)dr. 
We will now argue on how a can be identified with the 
surface tension of the light beam. In first place, we have 
computed numerically the parameter a for the different 
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high-power top-flat solitons, and we have found that it 
converges quickly to a fixed value a ~ 0.057 as soon as 
| /i | approaches the limiting value /i^. In such a limit, 
the central pressure goes to zero and the only important 
contribution to the integral defining a comes from the 
gradient term of Eq. ®. Thus a ~ J£° r\ V±^\ 2 dr. 

We can then get an analytical approximation for a by 
noting that for large R and r, ~ d^/dr, and the 

multiplicating r in the integrand can be approximated 
by R in the comparatively thin 'surface' region. Taking 
into account that for large r Eq. (pQ) yields A' (r) 2 ~ 
-7A(r) 4 /2 + (L4(r) 6 /3 + /M(r) 2 , and by substituting the 
limiting value \i 
the expression: 
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- , we get after some algebra 
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For instance, the choice 7 = J = 1 gives <r = 0.057, in 
complete agreement with our numerical value. By dif- 
ferentiating Eq. J5} and taking into account that a is 
constant, we get: 
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dR 
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Our numerical calculations show that the growth in R 
holds the thermodynamical equilibrium, (i.e.: ^| — 0), 
within a relative error which turns out to be as small as 
^ < 10" 4 § for all the range of "flat-top" eigenmodes 
considered. Therefore, we can set = in Eq. (fTDj). 

d(R 2 p c ) 
dR 



and we get 



2<T, or 
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(ii) 



which is the celebrated Young-Laplace (Y-L) equation [l3j 
for spherical liquid droplets being a the surface tension 
of the system. Therefore, despite the failure of the usual 
hydrodynamic approach, we have demonstrated that the 
pressure distribution in the inhomogeneous region has a 
deep physical interpretation. In fact, the integral of the 
pressure in the non-homogeneous surface region gives the 
surface tension (a), i.e. the inward force that compen- 
sates the outward positive inner force described by the 
pressure p c , in order to keep the droplet stationary. 

Instead of directly comparing Eq.([TT]) with the nu- 
merical simulation, we will equivalently test the inverted 
equation R — 2<r/p c , where p c can be expressed ana- 
lytically in terms either of the central amplitude A or 
of the propagation constant /i, by using Eq. ([5]) and 
Eq.([9]). This result provides the first analytical expres- 
sion for the radius R of the bidimensional top-flat soli- 
tons as a function e.g. of /i, and is compared with the 
numerical solutions of the stationary version of Eq.([Tj) in 
FigHJ As it can be appreciated in the figure, the agree- 
ment between the analytical formula and the numerical 
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FIG. 2: [Color online] Numerical simulation of the dripping 
behavior of light beams propagating through a channel waveg- 
uide highlighted with solid lines. As a comparison, the upper 
(lower) row corresponds to snapshots showing the evolution of 
the modulus of the quasi-gaussian ("flat-top") solutions. The 
size of the pictures in x and y is the interval [—250, 250]. The 
snapshots a)-d), b)-e), f) and c) correspond to propagation 
distances z = 0, 11000, 27000, 45000, respectively. 



computation is remarkable, and becomes complete when 
approaches \fioo\- This result confirms our theoretical 
framework and provides the first formal demonstration, 
by validation of the Y-L equation, of the liquid properties 
of the "flat-top" solutions in the (2 + 1) dimensional CQ 
model. Moreover, from Eq. ([TT]) it can be inferred that, 
as R — > 00, the value of p c vanishes, indicating that sur- 
face tension effects are not needed to balance the inner 
pressure, as it is the case of standard liquids described 
by the Y-L equation. 

Dripping of light droplets.- In classical fluid mechan- 
ics, the presence of surface tension effects can be appreci- 
ated in the dynamical phenomenon of capillarity [3, Q3] • 
Here, in order to study droplets formation and dripping 
in our system, we have introduced in our mathematical 
model an external "channel-type" linear optical waveg- 
uide V(x,y), superposed to the cubic-quintic nonlinear- 
ity. This waveguiding structure consists of three regions 
with indices m = 0.002 n 2 = 0.0028 and n 3 = 0.001, 
which correspond to the top region(ni), central channel 
and bottom zones (712) and rectangular regions flanking 
the central channel^), respectively, as it can be seen 
in the pictures of FigI2J In our simulations, we compare 
the evolution of two initial eigenstate beams with prop- 
agation constants /i//ioo = 0-01 (quasi-gaussian beam of 
Fig Hi) and fi/fi^ = 0.98 ("flat-top" beam of FigJUl), 
both located within the waveguide top region. Depend- 
ing on the channel size, and keeping 77,3 < ri2,^i, above 
a given value of An = n 2 — n\ in both cases a significant 
amount of beam power starts to flow from the initial 
eigenstate through the channel, as shown in (Figsl2j3 and 
[2^). It is noteworthy that the light stream inside the 
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channel does not suffer any unstabilization and remains 
connected to the initial source of light, i.e., the guide 
prevents the appearance of modulational instability. At 
the output of the channel, it can be seen in Figj2k that 
the low-power distribution spreads like a (coherent) gas 
in free expansion. However, the light flowing from the 
"flat-top" beam yields to the formation of a droplet, as 
it can be appreciated in Fig|2j and in more detail in the 
insets of FigOl 

As it can be seen in the insets of FigOl just before 
the release of the droplet the falling column of liquid 
light becomes narrower, just as in the well-known case 
of usual liquid streams falling from a tap[15]. On the 
other hand, all the steps in Figj2] are qualitatively very 
similar to those obtained both in real experiments and in 
simulations with liquids [I3]-[l7|. 
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FIG. 3: [Color online] Numerically calculated values of the 
surface tension along the droplet-formation process. Solid 
(dashed) line corresponds to numerical (analytical) calcula- 
tions of a. Insets: snapshots showing different stages of evo- 
lution of the light droplet for different values of z. Circles in 
the figure refer to upper insets. The spatial scales spanned 
are x £ [-50,50],?/ e [-20, -80]. 

Finally, we have also performed a quantitative test by 
comparing the numerical simulation with the generalized 
Y-L equation that describes the growth of an elliptic bub- 
ble of an ordinary liquid system, p c = cr(-^- + — ^- ), where 
Ri and Rn are the principal radii [13]. In Figj3[ we have 
plotted the quantity a = Pc/(j^ + j[jj) for the "dripping" 
simulation of Fig|2] at several propagation distances, as- 
suming z = where the light droplet first appears. We 
see that the numerical value of a oscillates around the 
same value that we have calculated analytically in Eq. ([9|) 
for the stationary solutions. This demonstrates that the 
droplets are formed close to stationary equilibrium, and 
they grow according to the generalized Y-L equation as 
in the case of an ordinary liquid. 

Conclusions.-. We have provided the first consistent 
computation of the pressure and surface tension of the 
soliton solutions appearing in the propagation of self- 



trapped laser beams described by the (2 + 1) dimensional 
CQ-NLSE, and we have shown both analytically and nu- 
merically that they satisfy the Young-Laplace equation 
that governs the equilibrium of usual liquid droplets. 
Subsequently, we have also demonstrated that the sys- 
tem dripping properties are governed by the same gen- 
eralized Y-L that applies to an ordinary liquid. These 
results reveal the deep connection between the physics 
of self-trapped laser beams and quantum liquids at zero 
temperature, opening the door for the quest of these new 
states of matter in the frame of current nonlinear optics 
experiments. 
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